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Particle tracking is a key to single-particle-level confocal microscopy observation of colloidal sus- 
pensions, emulsions, and granular matter. The conventional tracking method has not been able to 
provide accurate information on the size of individual particle. Here we propose a novel method to 
localise spherical particles of arbitrary relative sizes from either 2D or 3D (confocal) images either in 
dilute or crowded environment. Moreover this method allows us to estimate the size of each particle 
reliably. We use this method to analyse local bond orientational ordering in a supercooled poly- 
disperse colloidal suspension as well as the heterogeneous crystallisation induced by a substrate. 
For the former, we reveal non-trivial couplings of crystal-like bond orientational order and local 
icosahedral order with the spatial distribution of particle sizes: Crystal-like order tends to form in 
regions where very small particles are depleted and the slightly smaller size of the central particle 
stabilizes icosahedral order. For the latter, on the other hand, we found that very small particles 
are expelled from crystals and accumulated on the growth front of crystals. We emphasize that such 
information has not been accessible by conventional tracking methods. 



I. INTRODUCTION 

Experimental physicists often need to recognize ob- 
jects, count them, follow them and characterise them. 
Perrin [Ij had to count colloids by hand to establish 
the sedimentation-diffusion equilibrium. Nowadays com- 
puter vision algorithms are used routinely in the lab 
to track hundreds of thousands of objects as diverse as 
stars in a galaxy [2], tracers in a microfluidic device [3], 
pattern formation in polymer systems [H [5], dust in a 
plasma, bacteria |3 ^ or viruses in a living cell [8] . In 
all these cases, tracking is possible if the particles are ei- 
ther point-like and far apart, or several pixel wide and 
almost monodisperse in size (at least the smallest dimen- 
sion for anisotropic objects [6l |7]). To our knowledge, 
algorithm that allow tracking of polydisperse particles in 
crowded environments have not reached the soft matter 
community. 

Overall particle size distribution in colloidal suspen- 
sions and emulsions influences crystallisation [9HT3] , glass 
forming ability [9i [TOl [T2HT4] . sedimentation [I5l[l6] and 
emulsion stability [151 IIZ] among other physical phenom- 
ena. It can be characterized by various methods that 
rely on measurements done in well-controlled environ- 
ments [18-20 . However the local size distribution is not 
accessible experimentally in situ and has thus not been 
studied so far. 

Particle-level microscopy experiments usually access 
the coordinates of the particles via the algorithm pro- 
posed by Crocker and Crier [21]. The original noisy im- 
age is blurred by convolution with a Caussian kernel of 
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width cr to yield a soft peak per particle. Local intensity 
maxima within this blurred image give the coordinates 
of the particles with pixel resolution. Sub-pixel resolu- 
tion (0.1 0.3 pixels error) can be achieved by taking 
the centre of mass of a neighbourhood around the lo- 
cal maxima. The extension of this algorithm to localize 
particles in three-dimensional (3D) confocal microscopy 
images has been done in two ways: either tracking parti- 
cles in each confocal plane and reconstructing the results 
(2D-flavour) [22^, "23 , or full image analysis on three di- 
mensional pictures (3D-flavour) [24 . 

The choice of the width <j of the blurring kernel is 
critical: if it is too small, then the intensity profile is flat 
near the centre of a particle, leading to multiple and ill- 
localized maxima per particle; if it is too large, then the 
peaks of nearby particles overlap, leading to shifts in the 
detected positions [25l [26] , or even fusion of the particles 
(only one particle detected instead of two). If the colloids 
are fairly monodisperse one can argue (at least in the 3D- 
flavour) that there exists a range of possible width where 
the choice of a has almost no effect on the number of 
particles detected. Choosing a within this range gives 
confidence in the localisation results. 

However, we found that no such "good blur width" ex- 
ists in a sample of moderate (6 — 7%) polydispersity (see 
Fig. [ij^-h). The detection of smaller particles with small 
blurring is incompatible with the detection of the larger 
particles, and conversely. This unacceptable failure of the 
Crocker and Crier [2j algorithm, as well as the want of 
the particles' radii data, triggered our design of a novel lo- 
calisation algorithm that would be robust even for a sys- 
tem of any finite polydispersity, which is unavoidable and 
sometime desired in real experiments. Recently Kurita 
and Weeks [271 EH] have designed a sizing method using 
particle coordinates from confocal experiments. However 
their method do not work at the image processing level 
and relies on coordinates extracted via the Crocker and 
Crier pT algorithm. If these coordinates are wrong or 
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FIG. 1. Visualisation of the results of various tracking methods for the same portion of image, (a) Multiscale 3D tracking, (b) 
Reconstruction from multiscale 2D tracking, (c-h) Crocker and Grier method in 3D with blurring radius increasing from 2 px 
to 4.5 px by steps of 0.5 px. The circles on each picture are the result of 2D multiscale tracking of each XY slice of the 3D 
pictures. Sphere are displayed with radii determined by the tracking methods in (a) and (b), and equal to the blurring radius 
for (c)-(h). 



if some particles are missed (as shown in Fig. [T]), the 
output of their method could not be exact. 

The key notion to detect objects of unknown and pos- 
sibly diverse sizes in an image is the scale space [29 . A 
popular implementation for isotropic objects (or "blobs") 
is the Scale Invariant Feature Transform (sift) of Lowe 
[3Qj . It is often used to match between different images 
from complex objects consisting of many rigidly linked 
blobs {e.g., to create a large scale image from overlap- 
ping pictures) [3QH32] . To our knowledge, this method 
has never been used for the quantitative localisation and 
sizing of independent single-blob objects like spherical 
colloids, droplets in an emulsion, or crystal nuclei. 

Here we apply this new particle tracking method to 
study how the sizes of particles affect local structural or- 
dering in a supercooled colloidal suspension and the pro- 
cess of heterogeneous nucleation from a substrate. We 
reveal non-trivial local couplings between such orderings 
and the spatial distribution of particle sizes, which may 
provide crucial information for our understanding on how 
the polydispersity influences liquid dynamics and order- 
ing phenomena. 

The organization of our paper is as follows. In Sec- 
tion [nj we will describe our localisation method and its 
results on synthetic but more and more realistic data. 



In Section III of the paper, we will focus on the specific 
case of 3D confocal data. In Section [IV| we will apply 
our method to a crystallizing glass of polydisperse hard 
spheres observed by confocal microscopy and discuss the 
influence of size distribution on the ordering process. We 



conclude in Section M 



II. LOCALISATION AND SIZING METHOD 

In this section, we will start by recalling the principle of 
SIFT, then we will explain our method in the ideal case of 
an isolated binary ball, to add successively finite dilution 
and difference in brightness between the particles. 



A. Scale invariant feature transform 

The SIFT consists in convolving the original image / by 
Gaussian kernels g^^ of logarithmically increasing widths 
(Js to obtain a series of blurred images Gs 



V5>0, Gs=Ii^ga 



where ^ is the convolution operator and 



> 0, cr.s 



(1) 



(2) 



with n a fixed integer. Following Ref. [30 we use 
(Jo = 1.6 and n = 3. Bright objects in the original im- 
age appear as bright blobs in the blurred images, and 
the blobs fuse together as the kernel width increases (see 
Fig. [ij^-h). This can be seen as a series of low-pass fil- 
ters in the frequency domain. If we take the difference 
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FIG. 2. Application of sift in 2D. (a) Fluorescent droplets 
in a microfluidic device. The wall of the device is visible 
in the top-left corner, (b) Nucleation under phase contrast 
microscope. The blue circles are the result of our algorithm. 
For any tracking algorithm some detection failure near the 
edge of the image are unavoidable and indeed visible here 
since we show the whole pictures used for tracking. 



between consecutive blurred images we obtain a comb of 
band pass filtered versions of /: 

V8>1, DoGs = Gs-Gs-i (3) 

The difference of Gaussians {DoG) response function de- 
fined in this way depends on the position in space f and 
on the scale s. Bright objects in the original image are 
detected as local minima in DoG. With this procedure 
any feature with a radius as small as 2 px and as large 
as its distance to the edges of the image can be detected. 
Furthermore, the intensity of the response is optimal at a 
a that can be related to the size of the object (see below). 
Thus a local minima in both space and scale in the re- 
sponse function DoG indicates both localisation and size 
of an object, without any assumption on the target size. 

Because of the inherent polydispersity of many soft 
matter systems, the possible applications of SiFT-based 
localisation are countless. For example, droplets of an 
emulsion can be followed through a microfluidic device. 
Fig. [2^ shows the result of our version of the sift on a 
very polydispered and dense emulsion observed by wide- 
fleld fluorescence microscopy [33] . The sizes extracted are 
obviously correct except when the system cannot be con- 
sidered as 2D. Another possible application is nucleation 
rate measurement. When a phase transition proceeds 
via nucleation growth, the sizes of the nuclei are very di- 
verse, making their automatic counting difficult by other 
methods. As shown in Fig. the sift allows to count 
and to measure the size of nuclei during a liquid-liquid 
transition in a water-glycerol mixture [34 . 

However, to the extent of our knowledge, the sift has 
not been used in a physics context to obtain reliable mea- 
surements of particles numbers and sizes from well cal- 
ibrated images. We will address this reliability in the 
following. 

We will mainly cope with the three-dimensional exten- 
sion of the SIFT. As in the case of the Crocker and Grier 
[2T] algorithm, the SIFT can be extended in three dimen- 



sions {e.g., confocal microscopy images) in two ways: ei- 
ther by extracting 2D blobs independently on each slice 
and then reconstructing 3D objects (Fig. fib) or by work- 
ing directly in three dimensions [31, 32 (Fig. [l^). We 
found that the former method is prone to errors, missing 
about a tenth of the particles in our best implementation 
(Fig. [T]3). The 3D results presented in this paper are 
obtained solely by the later method. 

We stress that a volumetric implementation of sift im- 
plies a large amount of data (typically more than 1 Gb 
for a (256 px)^ picture) and thus requires careful mem- 
ory management. Our best implementation (C++) on a 
4 cores i7 computer takes less that 10 s to extract the po- 
sitions and scales of ~ 10^ particles in such a volumetric 
picture. A much slower implementation (Python+Scipy, 
single core) takes 1 s to deal with (1600 px)^ 2D images 
like Fig. We would expect real-time processing for 
2D images on GPU-enabled implementations. 



B. Sub-pixel and sub-scale resolution 

Here we assume perfect noiseless, distortion- free, im- 
ages. The objects to localise are thus (pixellised) balls of 
uniform intensity. To mimic the low resolution of exper- 
imental images in our test images (see Fig. [3|, we draw 
uniformly white balls on a 4 to 16 times larger image 
and then we reduce accordingly the resolution using area 
resampling. 

The DoG constructed above is defined on a (i + 1 di- 
mensional grid, with d the spatial dimensionality. One 
can interpolate it as a continuous function of the position 
r and the scale s and thus localise the minima with a pre- 
cision below the grid size, which corresponds to sub-pixel 
resolution on the position and sub-scale sizing [30 . We 
found that second order estimate of the spatial deriva- 
tives and first order estimates of the scale derivatives 
gave the best precision. 

The object-by-object optimal scale determination al- 
lows us to perform the spatial sub-pixel resolution step 
for each object on an image that is blurred just enough 
to have neither a fiat intensity profile nor an infiuence of 
the overlap with a nearby object's blob (an effect that 
plagues the Crocker and Grier [21 algorithm [25,. 26 J. 
We found that if for a given object the DoG is minimum 
at cTg, the best image to use is Gs-i- This leads to a spa- 
tial resolution below 0.1 pixels in the worst case (when 
particles are at hard-core contact), which is the same 
as the average precision claimed by Crocker and Grier 
[2T]. Moreover, when particle's surfaces are further than 
1 pixel, the error on the positions is less than 0.02 pixels 
(0.3% of the diameter) (see Fig. (3)3). 

C. Sizing at infinite dilution 

The analytical response of a binary ball to a Gaussian 
blur, at the centre of the ball is a function of dimension- 
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FIG. 3. Results from perfect images, (a-b) Sizing of an iso- 
lated sphere. Left of the vertical line our algorithm uses a 
doubled image, (c) Localisation error and (d) sizing error 
function of the distance between two particles. Oscillations 
are due to off-lattice centre position, (e) Size distribution 
extracted from digitized configuration of 4000 monodisperse 
hard spheres at 0.50 volume fraction (the vertical line indi- 
cates the input radius). The tail to the right is due to par- 
ticles on the edges of the image who have fewer neighbours 
and thus are more 'dilute', (d) and (e) also show the effect of 
finite dilution correction up to convergence. 



less ratio x = R/^/2a: G{x) (see Appendix for an exact 
expression). The DoG response function is the difference 
between the two such functions. However, the choice of 
the width of the two functions is not arbitrary: we make 
the difference between two consecutive blurred images, 
the image blurred by cTg+i and the image blurred by cTg. 



Therefore, in each of our discrete DoG images the value 
at the centre of the particle can be expressed as 

DoG{R, as, a) = DoG{xs,a) = G{xs/a) - G{xs), (4) 

with a = 2^/^. Given sub-scale refinement, this can be 
written as a continuous function of a: 

DoG{R, cr, a) = DoG{x, a) = G{x/a) - G{x). (5) 

Here it is clear that minimizing DoG{x, a) with respect 
to X yields a value x* that depends only on a. Exact 
calculation yields (see Appendix): 



din a 



1 



(6) 



where d is the spatial dimensionality. Practically one 
obtains <j*, the value of a that minimises DoG{R,(j,a), 
by a polynomial fit for discrete data aj . Eq. (|6| allows to 
translate the n-dependent a* to the parameter-free real 
radius of the particle, R. The error on R does depend on 
the number of subdivisions. We found that with n = 3 
the radius of an isolated pixelated ball can be indeed 
measured within 0.3% relative error with this method 
(see Fig. |3]d). 

The scale s = 0.5 corresponds (via Eq. ([2| and Eq. (|6|) 
to the smaller detectable radius Rmin ^ 3.5 px (ctq = 1.6, 
n = 3). In order to detect small objects, Lowe rec- 
ommends to double the size of the input image using 
linear interpolation prior to building the first level of the 
pyramid. This method fares relatively well in noiseless 
images (see Fig. [sJd) despite larger errors, but we found 
that it has to be used with care in confocal microscopy 
images due to noise and deconvolution artefacts. In ad- 
dition doubling an image size implies a 8-fold increase 
in memory consumption (reaching 60 GB in the case of 
a (512 px)^ original image). Our implementation allows 
this on 64 bit computers by relying on memory mapped 
files, nevertheless it is mostly impractical. All following 
tests results are obtained without relying on doubled im- 
ages. 



D. Edge and overlap removal 

The difference-of- Gaussian function has a strong re- 
sponse not only at the centre of bright objects but also 
along their edges. To eliminate these spurious detections, 
Lowe [30 suggests to construct the local Hessian matrix 
around each minimum of the DoG and then compare its 
eigenvalues to identify elongated objects. 

We found that this was often not enough, especially 
in crowded environments where an isolated void induces 
local minima of the DoG response with rather isotropic 
signatures on the edge of nearby particles. An other case 
not covered by the Hessian technique is the physically 
hierarchical structure of many soft matter systems, e.g., 
particles forming clusters. In this situation our algorithm 
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detects a blob for each particle and a much larger blob 
for the cluster. In both cases the DoG response of the 
spurious feature is smaller (less negative) than the one of 
the valid particles covering the same portion of space. 

Our method to remove the spurious feature is as fol- 
low: we looked for pairs of particles closer to each other 
than the larger of their radii (assuming infinite dilution). 
This means the centre of one of the particles is situated 
inside the other. In the physical systems we study, this 
cannot be correct due to excluded volume effects, thus 
we remove the feature of lesser DoG response. We name 
this method 'half overlap removal'. One may be tempted 
to implement a full overlap removal (no pair of parti- 
cle closer than the sum of their radii) to eliminate even 
more spurious detection. We found that the gain was 
extremely limited (0.01% of the detected features sup- 
pressed) and that imprecision in positioning or in sizing 
caused valid particles to be discarded. 

We note that the orientation of non-spherical particles 
or the deformation of soft objects can be monitored in 
principle by using the local Hessian matrix. 

E. Finite dilution 



where the matrix and the right hand side are computed 
given ({^4^}, {r^j}, cr*) and iteratively {Rj}^^\ with the 
upper parenthesised index indicating the iteration rank. 
The results of Eq. (|6| are good starting values for the 
radii. 

Using Eq. ([7|), the elements of the matrix simplify to 

In principle, Eq. (|9| is a x TV system of equa- 
tions. However the DoG functions and its derivatives 
are rapidly decaying functions, thus the matrices are ac- 
tually very sparse (about as many non-zero coefficients as 
particles in the first coordination shell), alleviating dra- 
matically the computational burden when using sparse 
system solvers. Fig. [3]i shows the result of such correc- 
tion for two identical particles, where Eq. ^ converges 
in a single iteration. In a many body case (see Fig. [3^) 
the convergence is reached in two iterations; however, 
the extremely low error of the dilute case is not totally 
recovered (~ 3% relative error rather than 0.3%). 



The Gaussian (and DoG) response of a particle de- 
cays rapidly away from its surface (see Eq. ( [21] ). We 
found that even at contact the position shift induced by 
a nearby particle was less than a tenth of a pixel (see 
Fig. [sj)). However when two particles are closer than 
a few times the blurring width, their infiuence on each 
other cannot be neglected: the minima of the DoG is ef- 
fectively shifted toward smaller scales, leading to smaller 
radii if one uses Eq. (|6| out of the dilute context (see 
Fig#). 

The response of N particles of radii {Ri} can be su- 
perimposed at any scale a and at any point of space, in 
particular we define DoGi as the response at the centre 
of particle i. 

DoGi{{Aj},{r,j},{Rj},a) = ^AjDoG{r = nj,R = Rj, 

3 

(7) 

where the function DoG{r^ i?, a) is the response of a par- 
ticle of radius R at distance r from its centre; r^j is the 
distances between particles i and j and {Ai} are the re- 
spective brightness of the particles. The sift algorithm 
yields {cr*} so that for all z, DoGi is minimum with re- 
spect to a at cr* , thus differentiating Eq. ^ with respect 
to a: 

Vi, ^^({A,},{r,,},{i?,},<T*)=0. (8) 

The system defined by Eq. ([8| is non-linear with re- 
spect to {Rj} but can be solved iteratively by Newton's 
method: 



F. Brightnesses 

To solve Eq. (|9|, one needs knowledge of the bright- 
nesses {Ai}. In a first approximation, they can be as- 
sumed equal to a constant, which allows to simplify them 
out. This is often a sensible approximation. Neverthe- 
less, the particles in an experimental image are not uni- 
formly bright due to synthesis imperfection (quantity of 
dye fixed by each particle) and photo bleaching. If one 
does not take into account the relative brightness of the 
particles, less bright particles will appear smaller. 

A better approximation is to measure during the SIFT 
process the value of the DoG response at the position and 
^^scale of each particle, i.e. DoGi{{Aj}^ {rij}^ {Rj}^a*). 
Given the (iterative) values of the one can solve 



^d^DoGi 



dRjda 



(fc+i) 



dDoGi 



, (9) 



Eq. ([7| to get an iterative value of {A^}*^^^. With respect 
to the brightnesses, Eq. ([7| is a linear system of equa- 
tions with N unknowns, thus directly solvable. It is also 
as sparse as Eq. 

To sum up, the coefficients {Ai} can be computed 
along with the radii in a joint iterative process: 

repeat 

until convergence 

In our tests, we found that for both cases with and 
without the brightness determination this algorithm con- 
verges quickly in one or two iterations. 
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III. APPLICATION TO 3D CONFOCAL 
MICROSCOPY IMAGES 

A. Effect of a point spread function 

Real images suffer from optical limitations, e.g.^ the 
point-spread function (psf) of the microscope. In par- 
ticular, spherical particles observed by confocal mi- 
croscopy appear elongated along Z. The influence of 
such anisotropic distortion on the Gaussian and DoG re- 
sponses is not trivial and most of the equations given in 
Appendix do not have analytical equivalents. In particu- 
lar, the minimum of the DoG response is found at larger 
values of R/a^ thus the naive use of the methods detailed 
in the previous section lead to overestimate particle sizes. 

In addition, the overlap of neighbouring particle im- 
ages is no longer negligible when the particles are aligned 
along Z. This leads to large imprecision in particle 
positions, especially in the case of anisotropic environ- 
ment (e.^., isolated pair of particles, interface, colloidal 
gel). In Fig. [4]i the radial distribution function of al- 
most monodisperse colloids with short range attraction 
displays a spurious shoulder before its first peak, indicat- 
ing that particles centres are often found closer than the 
sum of their hard core radii, as illustrated in Fig. ^p. 

We found that these issues can be better addressed by 
pre-processing the images (as detailed below) rather than 
post processing the sift output via analytical methods. 



Using Eq. (11) we obtain 



\^x[h]f 



(15) 



Here J^x indicates the Fourier transform along axis X. In 
point scanning confocal imaging, the PSF has negligible 
lobes along X and Y {X only for line scanning), thus we 
have 



\^z[ 



l-F; 



(16) 



For example an Hemiitian kernel (real valued spectrum) 
is 



^z[h] 



\:fzW 



(17) 



Fig. ^p-c shows the particles localized from the original 
image and from the deconvolved image. Deconvolution 
mends both size overestimation and imprecision in z co- 
ordinate. This also translates in the radial distribution 
function (Fig. [Iji), where the spurious shoulder before 
the first peak disappears. 



Deconvolution 



C. From optical to hard-core radius 



The image y acquired by a microscope can be expressed 



as 



y 



X h -\- e. 



(11) 



where x is the perfect image, h is the PSF of the micro- 
scope and e is the noise independent of both x and h. 
The process of estimating x from y and some theoreti- 
cal or measured expression of h is called deconvolution. 
Deconvolution in the presence of noise is a difficult prob- 
lem [35 . Hopefully here we do not need to reconstruct 
the original image, but only our first Gaussian blurred 
version of it, i.e.. Go = ^^^ao starting from yo = yi^gao- 
Indeed, after a reasonable amount of blur in three di- 
mensions, the noise can be neglected and we thus obtain: 



^0 - Go ^ h, 



or in Fourier space 



(12) 



(13) 



Once T[h] is known the deconvolution reduces to a 
simple division in Fourier space. Let us measure T[h] in 
an isotropic system where we can write 



(14) 



The method exposed above relies on the uniformity of 
the brightness within each particle to extract optical radii 
that is a good approximation to the physical (hard core) 
radii. However a colloidal particle often shows a smooth 
intensity profile in itself under the microscope, less bright 
at the edge than at the centre of the particle. This can 
be due to inhomogeneous fixation of the dye during the 
colloid synthesis, but the bottom line is the in-plane PSF 
not corrected for in our above deconvolution method. 

Under such conditions, the particles are detected 
smaller than their expected sizes. The real to optical 
size ratio has to be set using our knowledge of the sam- 
ple. For example one may use the position of the first 
peak of the g{r) to measure the average hard-core diam- 
eter. In general the real to optical size ratio may depend 
on the size of the particles and may evolve in time (i.e., 
because of photobleaching). We present below a detailed 
analysis of these issues for a case system. 

We also note that darker edges lead to smaller influ- 
ence of a particle on the DoG response of its neighbours. 
A smaller optical radius makes the particles optically fur- 
ther from each other. Accordingly, we found that finite 
dilution corrections were less important in real images 
than in our synthetic test images. For the same reason, 
the scaling factor must be measured and applied after — 
not before — finite dilution corrections. 



FIG. 4. Deconvolution. Detail of the same YZ slice of (a) original confocal image, (b) previous blurred by ctq = 1.6, (c) 
previous deconvolved by measured kernel. Circles indicate the tracked particles position and size when using whether (b) or (c) 
as first Gaussian layer Gq. All three centres are in the slice ±0.5 px. (d) Radial distribution function of almost monodisperse 
sticky spheres localised (blue) without and (red dash) with deconvolution. 



IV. REVEALING POLYDISPERSITY EFFECTS 
ON STRUCTURAL ORDERINGS IN A GLASSY 
COLLOIDAL HARD SPHERE LIQUID 

Here we apply our tracking method to access couplings 
between the spatial particle size distribution and local 
structural ordering in a supercooled colloidal liquid with 
size polydispersity. 



A. Experimental 

We used PMMA (poly (methyl methacrylate)) colloids 
sterically stabilized with methacryloxypropyl terminated 
PDMS (poly (dimethyl siloxane)) and fluorescently la- 
belled with rhodamine isothiocyanate chemically bonded 
to the PMMA. The colloids were suspended in a sol- 
vent mixture of cis-decalin and cyclohexyl-bromide for 
both optical index and density matching. To screen 
any (weak) electrostatic interactions, we dissolved tetra- 
butylammonium bromide salt, to a concentration of 
300 nmol/L [36 . The estimated Debye screening length 
is 13 nm, well below the length scale of the colloids 
that can be considered as hard spheres. The 8 bit 
graylevel data was collected on a Leica SP5 confocal mi- 
croscope, using 532 nm laser excitation and voxel size of 
(283 nm)2 x 293 nm. 

To realize a more precise density matching, the tem- 
perature was controlled on both stage and objective lens. 
This setup allows us to alter the buoyancy simply by a 
temperature change, thus to set the effective gravity up- 
ward, downward or almost null. After careful shear melt- 
ing, the sample was filled into a 100 iitm x 1 mm capillary 
(Vitrocom) and set on the microscope stage. We then 
spent a few days to find the temperature corresponding 
to exact density matching (within 0.1 K). This waiting 
time was enough for crystallites to form on both top and 
and bottom walls. Then, we heat up our sample by a 
few degrees compared to the density-matching tempera- 
ture in order to make the colloids heavier than the sol- 



vent. At the top of the sample (close to the objective 
lens and thus allowing much clear imaging) the volume 
fraction drops and all crystallites melt. Finally, we set 
back the thermostat to the density-matching tempera- 
ture, allowing the top of the sample to slowly return to 
its supercooled state. We could then observe the hetero- 
geneous nucleation from the beginning. We tracked the 
sample from top to bottom and we are thus sure that 
the crystallite actually form at the wall and do not come 
from the rest of the sample. Namely, crystallisation in 
this sample is caused only by heterogeneous nucleation 
on the substrates, presumably due to the wall-induced 
enhancement of crystal-like bond orientational ordering 



B. Global size distribution 

We localise the particles using our method with a pre- 
blur of only <Jo = 1.0 (see Section II A) to be able to 
detect the smaller particles {Rmi 



^ 2 px) without ex- 
pensive oversampling. We checked that choosing higher 
(Jo truncates the size distribution but has no other effect 
on its shape. After removing half-overlapping features, 
we applied a single iteration of finite dilution corrections 
for both intensities and radii to obtain the optical radius 
of each particle. 

We then estimate the hard-core diameter function of 
the optical radius R by locating the first peak of the par- 
tial radial distribution function gnir) of the particles hav- 
ing Ri close to R. As shown in Fig.|5^, the real to optical 
size ratio is around 1.5 and is rather constant respective 
to i?, thus a single overall scaling factor is enough. It 
can be determined by a single partial gR{r) with R near 
the peak of the size distribution, or by locating the peak 



of ^(r), with 



ij/{Ri ^ Rj). We found that the 



real-to-optical size ratio was increasing with time due to 
photobleaching (see Fig. ^p). We fit this increase by a 
linear relation and applied the resulting time-dependent 
ratio to obtain the real size of each particles. 
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The consistency of our method can be checked by 
constructing the radial distribution function g{r). In 
monodisperse hard spheres, the g{r) should have a sharp 
first peak at r = 2R corresponding to hard core contact. 
Polydispersity implies hard core contacts at various r and 
thus broadens the peak. One can recover a sharp peak by 
constructing g{r) (this time r is scaled by the real radii, 
not the optical ones). In Fig. [sji we successfully used the 
sizes measured by our method to rectify the first peak. 

Fig. ^ shows the size distribution only 140 dry col- 
loids measured on high resolution scanning electron mi- 
croscopy (sem) images. We also show the size distribu- 
tion obtained by our in situ measurements (~ 1.7 x 10^ 
instantaneous sizing). The main peak of the later com- 
pares well with the former once a solvent swelling of 25% 
in radius is taken into account. However the small sam- 
pling of the SEM measurements completely misses the 
tail toward small sizes featuring two low peaks that prob- 
ably correspond to secondary and ternary nucleation dur- 
ing the synthesis of the colloids [38l [39] . 

From our data we compute an overall volume fraction 
of 0.60, almost constant during the experiment. The 
polydispersity estimated from SEM data is 6.2%. This is 
coherent with a fit of the main peak of our data (6.9%), 
however the polydispersity of the whole distribution in- 
cluding the tail toward small sizes is 14.8%. We will see 
below how such a complex size distribution affects the 
physical behaviour of the system. 



C. Polydispersity and structural heterogeneity in a 
supercooled colloidal liquid 

Hard sphere supercooled liquids and glasses contain 
both medium ranged crystal-like bond orientational or- 
dering (mrco) and local icosahedral ordering [40 . The 
later acts as a frustration against the expansion of the 
former and thus against crystallisation. Both kind of 
structures can be detected using bond orientational order 
(boo) introduced by Steinhardt et al. [4l (see Fig. [6|. 
As we have described elsewhere [40^, crystal-like bond 
ordering is well described by the scalar bond order pa- 
rameter Qe, and the icosahedral bond ordering by wq. 
FCC or HCP crystals under thermal vibrations typically 
have Qq > 0.4 [42]. We distinguish MRCO from liquid 
structures by a threshold at Qe = 0.25. The perfect 13 
particles icosahedron is the minimum of wq^ at negative 
values. We consider that a neighbourhood is icosahedral 
when wq < —0.023. 

In Fig. [7| we compare the overall size distribution to 
the size distribution within each remarkable structure. 
The roles of the particle at the centre of an icosahedron 
and of the particles surrounding it are asymmetric. The 
size distribution of the particles at the surface of icosa- 
hedra is almost identical to the overall size distribution, 
however the size distribution of the particles at the cen- 
ter of icosahedra features a second peak at radii about 
80% smaller than the main peak. A centre to surface size 
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FIG. 5. In situ sizing of colloids in a glass, (a) Position of the 
first peak of the partial gnir) function of the optical radius R. 
Solid line corresponds to a real-to-optical size ratio of 1.5. (b) 
Time dependence of the real-to-optical size ratio. Solid line is 
a linear fit. (c) Size distribution estimated by our algorithm 
(dashed line). Comparison with the estimation from SEM of 
only 140 dry particles (steps) is possible once 25% of swelling 
is taken into account (full line), (d) First peak of the radial 
distribution function with (full line) and without (dashed) the 
individual sizes data. 



ratio near 0.8 seems to stabilize icosahedral order. This 
is consistent with a recent simulation study by Shimono 
and Onodera [43] . 

We note that the tail of small particles characteristic of 
our overall size distribution is less pronounced in MRCO 
and almost disappears in well- formed crystals (Fig. [7^) . 
The width of the main peak of the distribution does 
not change significantly upon crystal-like ordering. This 
means that a small amount of polydispersity is acceptable 
inside MRCO, coherently with the results for the (defec- 
tive) crystalline lattice by Fasolo and Sollich [llj. How- 
ever MRCO cannot appear where markedly smaller parti- 
cles are present. We also checked the difference between 
FCC {w4 < 0) and HCP {w4 > 0) to find that they have 
both exactly the same size distribution (not shown). This 
suggests that the size distribution controls the formation 
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FIG. 6. Structure visualisation at t = (left) and t = 40 h 
(right). Small particles {R < 1.5 ^im) are shown in green, 
crystal-like ordered particles {Qq > 0.25) in red and icosahe- 
dral particles and their neighbours in purple. Other particles 
are not shown. 
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FIG. 7. Size distribution within local structures in a super- 
cooled hard sphere colloidal liquid, (a) With increasing local 
crystalline bond orientational order: all particles, MRCO par- 
ticles (0.4 > Qe > 0.25) and crystalline particles {Qq > 0.4). 
Crystalline ordered regions are less polydisperse. (b) Particles 
at the centre or on the surface of icosahedra {wq < —0.023). 
A small central particle with large surface particles tend to 
stabilize icosahedra. 



of MRCO but has no influence on its symmetry (proba- 
bly, the crystal polymorph). We stress that this coupling 
of MRCO and icosahedral ordering with the spatial dis- 
tribution of particle sizes does not imply fractionation, 
but rather suggests that a rather homogeneous size dis- 
tribution tends to help or stabilize MRCO and a slightly 
smaller size of a particle also stabilizes the formation of 
icosahedral order around it. We will see below that the 




Z/ (2i?peak) Z/ (2i?peak) 



FIG. 8. Instantaneous density profiles at t = (left) and 
t = 40 h (right). Top: large particles {R > 1.5 ^tm) with- 
out smoothing (black line) and with a Gaussian smooth- 
ing of width 2i?peak (gray area). Bottom: small particles 
{R < 1.5 ^im) with the same color code. For all figures, the 
red dashed curve is the smoothed density of ordered particles 
(Qe > 0.25) irrespective of their size. 



fractionation appends later in the crystallisation process. 



D. Polydispersity and heterogeneous crystallisation 
in a supercooled colloidal liquid 

As explained above, we triggered heterogeneous crys- 
tal nucleation on the top wall of our container. At each 
time step we can construct a z-dependent density profile 
for any species of interest. In Fig. [S] we show at two time 
steps the density profiles of large particles {R > 1.5 M-m), 
small particles {R < 1.5 iitm) and particles with medium 
to high crystalline order {Qq > 0.25). Icosahedral par- 
ticles (not shown) are an order of magnitude fewer than 
small particles (see Fig. [6| and are basically irrelevant. 
Near the walls, the density profiles show oscillations char- 
acteristic of layering (see, e.g.^ Refs. SID^ whose 
wavelength corresponds to the peak of the size distri- 
bution (the 'majority species'). We note that the small 
particles also shows layering on the same wavelength, al- 
though with a smaller amplitude. 

Once the short oscillations removed via a Gaussian 
smoothing of width 2i?peak5 one can see that the den- 
sity of large particles is constant from top to bottom of 
the sample (the slight dip on the edges is due to a small 
misalignment between our z axis and the normal of the 
wall). It confirms the accuracy of our density matching 
but also stresses that crystallisation takes place without 
an increase in density. The smoothed density profile of 
small particles starts almost flat. However the small par- 
ticles are expelled from the forming crystal, and as the 
crystalline front propagates it pushes in front of itself a 
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FIG. 9. Heterogeneous nucleation. Successive computer re- 
construction of the first three layers from the wall, every 12 h, 
seen from the fluid side. Particles sizes are according to our 
method. Small {R < 1.5 ^im) particles are shown in green. 
Other particles are coloured according to their crystal-like 
order. Note how the small particles are expelled from the 
forming crystallites to concentrate in the grain boundaries. 



growing ridge of small particles, very apparent after 40 h 
in both Fig. [8] and Fig. [6| 



V. CONCLUSION 

We have presented a new method to extract the coor- 
dinates of particles of very diverse sizes even when they 
are very close (at contact) to each other. This opens new 
experimental possibilities in soft matter systems where 
the particles are polydisperse {e.g., colloids, emulsions, 
and granular particles) or even changing size {e.g., phase 
separation dynamics). Our method also reliably extracts 
the size of individual particle, which allows a in situ anal- 
ysis of the consequences of naturally occurring size dis- 
tributions (beyond monodisperse, bidisperse or Gaussian 
distribution). 

We showed that in the case of colloidal hard sphere, a 
size distribution with a main peak and a tail toward small 
sizes induces different behaviours than what would be 
expected from a Gaussian size distribution. Our method 
allowed us to measure this distribution properly and to 
look into the interactions between size distribution, local 
ordering and heterogeneous crystal nucleation in a glass. 
We found that the presence of a long tail toward small 
sizes in the size distribution, was leading to highly poly- 
crystalline materials that may be of interest for engineer- 
ing purposes. We hope that this novel particle tracking 
method with a capability of particle-size determination 
will be applied to a wide field in soft matter physics. 



We confirmed this view by looking at the behaviour of 
small particles and crystalline order in a constant plain 
(Fig. [9] and Supplementary Movie 1). Before crystalli- 
sation, small particles are scattered randomly and MRCO 
develop where few very small particles are present. Then, 
crystallisation proceeds from the MRCO, pushing away 
the small particles that accumulate in the grain bound- 
aries. As a consequence, heterogeneous crystallisation do 
not proceed via a spatially homogeneous layer-by- layer 
growth, as observed at low polydispersity [45], but by 
many nucleation-growth events reflecting the initial fluc- 
tuations of the MRCO. This dynamics is reminiscent of 
the homogeneous nucleation case [46H48] but enhanced 
by the wetting of the wall by crystalline bond orienta- 
tional order [37 . 

To sum up, the small particles affect the growth of the 
crystal by coupling ordering to diffusion. Without the 
small particles, moving grain boundaries implies mov- 
ing particles only locally and thus propagation is easy. 
With very small particles, moving grain boundaries im- 
plies moving them by the same amount by diffusion. In 
the same way, the crystallisation front has to push the 
small particles in front of it, effectively increasing the 
polydispersity in the melt and thus slowing its propa- 
gation. Indeed, weeks after preparation our samples was 
still only partially crystallised and highly polycrystalline. 
We note that this coupling is due to the conserved nature 
of the species 'small particles'. Non-conserved species 
defined by local structural ordering such as icosahedra 
cannot produce such a coupling. 
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APPENDIX: GAUSSIAN RESPONSE OF A 
BINARY BALL 

A binary ball Br of radius R is convolved by a Gaus- 
sian kernel of width a. The response at a distance r 
(along z-axis without loss of generality) is 

G{r,R,(j) =-——^j7^ I e 2.^2 dxdydz 

(18) 
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or in spherical coordinates (p, 0, (/)) after integration along The difference of Gaussians response is 



27r<j^ Jo 
1 



p^e-^dp I sin l9e^"°^^dl9 

(19) 



DoG{r, R, cr, a) = G{r, R, aa) - G{r, R, a) (28) 



rcr v27r 



pe 



io 



(20) 



this can be integrated using the error function erf 



and has the partial derivative relative to the scale 



G(r, J?, ct) = - (5(r, R, a) + g{-r, R, a)) 



with 



g{r,R,a) = erf 



R^r 



2 a _ (fi±r)! 



(21) 



(22) 



dDoG, ^ , dG, ^ ^ dG, ^ , 

{r, R, cr, ol) = a-^{r, R, aa) - (r, R, a) (29) 



At the centre of the ball, Eq. pT]) reduces to 
R 



G{r = 0,R,(T) =erf 



aV2 



2 R _Rl_ 

e 2^x2 

TT a 



(23) 



In addition, we compute the following useful partial 
derivatives 



dg_ 
da 

dadR 
dG 



{r,R,a) = 



{R^ + rR + a^) _(r± 



{r,R,a) = -R 



da 
dadR 



{r = 0,R,a) 



ir = 0,R,a) 



^ '-^= — e 2^2 

2R^ _Bl_ 



TT (7^ 

2 .i?2_3^2 

R a ^ 



(24) 
(25) 
(26) 
(27) 



When the difference of Gaussians is minimum at the cen- 



tre of the ball, we combine Eq. (26) and Eq. (29) to get 



3 In a 



R 



Ty/2 V 1 - a 



(30) 



One can show that the factor 3 in the radical is the 



dimensionality of the space, thus Eq. (30) can be gener- 
alised to other dimensions. 
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